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Abstract 



This paper is concerned with exact real solving of well-constrained, bivariate polynomial systems. 
The main problem is to isolate all common real roots in rational rectangles, and to determine their 
^ , intersection multiplicities. We present three algorithms and analyze their asymptotic bit complexity, 

f^^ ■ obtaining a bound of Ob{N^'^) for the purely projection-based method, and Ob{N^^) for two subresult- 

ant-based methods; this notation ignores polylogarithmic factors, where A'^ bounds the degree and the 
bitsize of the polynomials. The previous record bound was Ob{N^'^). 

Our main tool is signed subresultant sequences. We exploit recent advances on the complexity of 
^f^ \ univariate root isolation, and extend them to sign evaluation of bivariate polynomials over two algebraic 

C^ . numbers, and real root counting for polynomials over an extension field. Our algorithms apply to the 

O^ ' problem of simultaneous inequalities; they also compute the topology of real plane algebraic curves in 

Ob{N^^), whereas the previous bound was C's(A'^^*). 

All algorithms have been implemented in maple, in conjunction with numeric filtering. We compare 
them against fgb/rs, system solvers from SYNAPS, and maple libraries insulate and TOP, which com- 
J\. , pute curve topology. Our software is among the most robust, and its runtimes are comparable, or within 

?— j ' a small constant factor, with respect to the C/C-l— I- libraries. 

cd ■ 
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1 Introduction 

The problem of well-constrained polynomial system solving is fundamental. However, most of the algorithms 
treat the general case or consider solutions over an algebraically closed field. We focus on real solving 
of bivariate polynomials in order to provide precise complexity bounds and study different algorithms in 
practice. We expect to obtain faster algorithms than in the general case. This is important in several 
applications ranging from nonlinear computational geometry to real quantifier elimination. We suppose 
relatively prime polynomials for simplicity, but this hypothesis is not restrictive. A question of independent 
interest, which we tackle, is to compute the topology of a real plane algebraic curve. 

Our algorithms isolate all common real roots inside non-overlapping rational rectangles, output them 
as pairs of algebraic numbers, and determine the intersection multiplicity per root. Algebraic numbers are 
represented by an isolating interval and a square-free polynomial. 

In this paper, Ob means bit complexity and Ob means that we are Jgnoring polylogarithmic factors. 
We derive a bound of OBiN^"^), whereas the previous record bound was Ob{N^^) [16|, see also [3|, derived 
from the closely related problem of computing the topology of algebraic curves, where A'' bounds the degree 
and the bitsize of the input polynomials. This approach depends on Thom's encoding. We choose the 
isolating interval representation, since it is more intuitive, and is used in applications. In [16| . it is stated 
that "isolating intervals provide worst [sic] bounds". It is widely believed that isolating intervals do not 
produce good theoretical results. Our work suggests that isolating intervals should be re-evaluated. 

Our main tool is signed subresultant sequences (closely related to Sturm-Habicht sequences), extended to 
several variables by binary segmentation. We exploit the recent advances on univariate root isolation, which 
reduced complexity by one to three orders of magnitude to Ob{N^) [11|, |12|, [15|. This brought complexity 
closer to ©^(Af'^), which is achieved by numerical methods |26| . 

In [l9[, 2x2 systems are solved and the multiplicities computed under the assumption that a generic 
shear has been obtained, based on [32|. In [36[, 2x2 systems of bounded degree were studied, obtained as 



projections of the arrangement of 3D quadrics. This algorithm is a precursor of ours, see also IJ], except that 
matching and multiplicity computation was simpler. In [23|, a subdivision algorithm is proposed, exploiting 
the properties of the Bernstein basis, with unknown bit complexity, and arithmetic complexity based on 
the characteristics of the graphs of the polynomials. For other approaches based on multivariate Sturm 



sequences the reader may refer to e.g. [22|, |27 




For determining the topology of a real algebraic plane curve, the best bound is OBiN^^) [3|, ll6|. In |37| 
thr ee p rojections are used; this is implemented in INSULATE, with which we make several comparisons. Work 

13| offers an efhcient implementation of resultant-based methods, whereas Grobner bases are employed 
To the best of our knowledge, the only result for topology determination using isolating intervals is 

where a Ob{N^^) bound is proved. 

We establish a bound of Ob{N^'^) using the isolating interval representation. It seems that the complexity 
in [l6| could be improved to OBiN^'^) using fast multiplication algorithms, fast algorithms for computations 
of signed subresultant sequences and improved bounds for the bitsize of the integers appearing in computa- 
tions. To put our bounds into perspective, the input size is in Ob(A^'^), and the total bitsize of all output 
isolation points for univariate solving is in Ob{N'^), and this is tight. Notice that lower bounds in real 
algebraic geometry refer almost exclusively to arithmetic complexity [5| . 

The main contributions of this paper are the following: Using the aggregate separation bound, we im- 
prove the complexity for computing the sign of a polynomial evaluated over all real roots of another (Lemma 
[7]). We establish a complexity bound for bivariate sign evaluation (Theorem [H]) . which helps us derive 
bounds for root counting in an extension field (Lemma I2ip and for the problem of simultaneous inequalities 
(Corollary [M]). We study the complexity of bivariate polynomial real solving, using three projection- based 
algorithms: a straightforward grid method (Theorem [TS)) . a specialized RUR (Rational Univariate Represen- 
tation) approach (Theorem [T9| , and an improvement of the latter using fast GCD (Theorem [20| . Our best 
bound is OBiN^^); within this bound, we also compute the root multiplicities. Computing the topology of 
a real plane algebraic curve is in Ob{N^^) (Theorem [25]) . 

We implemented in MAPLE a package for computations with real algebraic numbers and for implementing 
our algorithms. It is easy to use and integrates seminumerical filtering to speed up computation when the 
roots are well-separated. It guarantees exactness and completeness of results; moreover, the runtimes are 



quite encouraging. We illustrate it by experiments against well-established C/C++ libraries fgb/rs and 
SYNAPS. We also examine MAPLE libraries INSULATE and TOP, which compute curve topology. Our software 
is among the most robust; its runtime is within a small constant factor with respect to the fastest C/C++ 
library. 

The next section presents basic results concerning real solving and operations on univariate polynomials. 
We extend the discussion to several variables, and focus on bivariatc polynomials. The algorithms for 
bivariate solving and their analyses appear in Section |4j followed by applications to real-root counting, 
simultaneous inequalities and the topology of curves. Our implementation and experiments appear in Section 

A preliminary version of our results appeared in [9| . 

2 Univariate polynomials 

For / e Z[yi,..., yk,x], dg(/) denotes its total degree, while (ig^{f) denotes its degree w.r.t. x. C{f) 
bounds the bitsize of the coefficients of / (including a bit for the sign). We assume Ig (dg(/)) = 0{C (/)). 
For a e Q, £ (a) is the maximum bitsize of numerator and denominator. Let M (t) denote the bit complexity 
of multiplying two integers of size t, and M (d, r) the complexity of multiplying two univariate polynomials 
of degrees < d and coefficient bitsize < r. Using FFT, M (r) = (5b(t) and M {d, t) ^ OBidr). 

Let /,g G I.[x], dg(/) = p > q = dg{g) and C{f),C{g) < t. Wc use rem {f,g) and quo(/,.9) for the 
Euclidean remainder and quotient, respectively. The signed polynomial remainder sequence of /, g is i?o = /, 
i?i = g, i?2 = — rem (/, .g), . . . , Rk ^ — rem {Rk-2,Rk-i), where rem {Rk~i,Rk) = 0. The quotient sequence 
contains Qi = quo {Ri, Ri+i), i = . . . fc — 1, and the quotient boot is (Qo, • • • , Qk~i, Rk)- 

We consider signed subresultant sequences [SLwhich contain polynomials similar to the polynomials in 
the signed polynomial remainder sequence; see [35| for a unified approach to subresultants. They achieve 
better bounds on the coefficient bitsize and have good specialization properties. In our implementation we 
use Sturm-Habicht (or Sylvester-Habicht) sequences, see e.g. [sj, [iSl, |20[ . By SR(/, g) we denote the signed 
subresultant sequence, by sr(/, g) the sequence of the principal subresultant coefficients, by SRQ(/, 5) the 
corresponding quotient boot. By SRj (/,(;), or simply SRj if the corresponding polynomials can be easily 
deduced from the context we denote an element of the sequence; similarly for sr^ and SRQ . Finally, by 
SR(/, (7; a) we denote the evaluated sequence over a S Q. If the polynomials are multivariate, then these 
sequences arc considered w.r.t. x, except if explicitly stated otherwise. 

Proposition 1. j^d. \23iJ Assuming p > q, SR(/, g) is computed in Ob{p^Qt) and £(SRj(/, g)) = 0{pt). 
For any f,g, their quotient boot, any polynomial in SR(/, 5), their resultant, and their gcd are computed in 
Osipgr). 



The following proposition is a slightly modified version of the one that appeared in |2d. |28|. 



Proposition 2. Let p> q. We can compute SR(/, g; a), where a G Q U {±00} and C (a) ~ a, in OB{pqT + 
q^a+p^a). ///(a) is known, then the bound becomes OBipqT + q^a). 

Proof. Let SRg+i = / and SRg ~ g. For the moment we forget SR^+i. We may assume that SRg-i is 
computed, since the cost of computing one element of SR(/, g) is the same as that of computing SRQ(/, g) 
(Pr. [1]) and we consider the cost of evaluating the sequence SR(5, SRg_i) on a. 

We follow Licktcig and Roy [20| . For two polynomials A, B of degree bounded by D and bitsize bounded 
by L, we can compute SR(A, _B;a), where £ (a) < L, in Ob(M {D,L)). In our case D = 0{q) and L = 
0{pT + qa), thus the total costs is Osipqr + q^cr). 

It remains to compute the evaluation SRg+i(a) = /(a). This can be done using Horncrs' scheme in 
C'B(pmax{T,pcr}). Thus, the whole procedure has complexity 

Osipqr + q^a +pmax{T,pCT}), 
where the term pr is dominated by pqr. D 



When q > p, SR(/, g) is f,g, — /, —{g mod (— /)) • ■ • , thus SR(/, 5; a) starts with a sign variation irre- 
spective of sign((/(a)). If only the sign variations are needed, there is no need to evaluate g, so Proposition [2] 
yields ObJ'PQt + p'^o')- Let L denote a list of real numbers. VAR(L) denotes the number of (possibly modified, 
see e.g. |J, lla, [l^) sign variations. 

Corollary 3. For any f,g, VAR(SR(/, g; a)) is computed in OsipqT + m.m.{p,q}^a), provided sign(/(a)) is 
known. 

We choose to represent a real algebraic number a £ KaZg by the isolating interval representation. It 
includes a square-free polynomial which vanishes on a and a (rational) interval containing a and no other 
root. By fred we denote the square-free part of /. 



Proposition 4. fli . \l2i . \la l Let f G Z[x] have degree p and bitsize Tf. We compute the isolating interval 
representation of its real roots and their multiplicities in Ob{p^ +P^''"f)- The endpoints of the isolating 
intervals have bitsize 0{p^ + pTf) and C (fred) = 0{p + Tf). 

Notice that after real root isolation, the sign of the square-free part fred over the interval's endpoints, 
say [a, b] is known; moreover, fred{^)fred{^) < 0. The following proposition takes advantage of this fact and 
is a refined version of similar proposition in e.g. [3|, Il5| . 

Corollary 5. Given a real algebraic number a = (/, [a,b]), where C{a) = ^(b) = 0{pTf), and g G Z[.t], 
such that Ag{g) = q and bitsizeg = Tg, we compute sign(g{a)) in bit complexity 03(^5 maxJTj, Tg} + 
pmin{p,q}'^Tf). 

Proof. Assume that a is not a common root of / and g in [a, b], then it is known that 
sign(g(a)) = [VAR(SR(/,g;a)) - VAR(SR(/,g; b))] sign(/'(a)). 



Actually the previous relation holds in a more general context, when / dominates g, see 38| for details. 
Notice that sign(/'(a)) = sign(/(b)) — sign(/(b)), which is known from the real root isolation process. 

The complexity of the operation is dominated by the computation of VAR(SR(/, g; a)) and VAR(SR(/, g; b)), 
i.e. we compute SRQ and evaluate it on a and b. 

As explained above, there is no need to evaluate the polynomial of the largest degree, i.e. the first (and 
the second ii p < q) of SR(/, 9) over a and b. The complexity is that of Corollary [3J i.e. OB(p'3'niax{T/,Tg}-t- 
rBin{p, q}^pTf). Thus the operation costs two times the complexity of the evaluation of the sequence over 
the endpoints of the isolating interval. 

If a is a common root of / and g, or if / and g are not relative prime, then their gcd, which is the last 
non-zero polynomial in SR(/, g) is not a constant. Hence, we evaluate SR on a and b, we check if the last 
polynomial is not a constant and if it changes sign on a and b. If this is the case, then sign{g{a)) — 0. 
Otherwise we proceed as above. D 

Proposition |4] expresses the state-of-the-art in univariate root isolation. It relies on fast computation of 
polynomial sequences and the Davenport-Mahler-Mignotte bound, see [8| for the first version of this bound. 
The following lemma, a direct consequence of Davenport-Mahler-Mignotte bound, is crucial. 

Lemma 6 (Aggregate separation). Given f G 1j[x], the sum of the bitsize of all isolating points of the real 
roots of f is 0{p^ +pTf). 

Proof. Let there be r < p real roots. The isolating point, computed by a real root isolation subdivision 
algorithm lll.ll2l Il5|. between two consecutive real roots, say aj and a^+i, is of magnitude at most 2\'^j ~ 



aj+i\ := ^Aj. Thus their product is ^ 0^=1 ^j- Using the Davenport-Mahler-Mignotte bound, the product 
is bounded from below, that is J| Aj > 2~'^^p +P^f) _ Taking logarithms, we conclude the proof. D 

We present a new complexity bound on evaluating the sign of a polynomial g{x) over a set of algebraic 
numbers, which have the same defining polynomial, namely over all real roots of f{x). It suffices to evaluate 
SR(/, g) over all the isolating endpoints of /. The obvious technique, e.g. [15[, see also [3J, |31|, is to apply 
Corollary [5] r times, where r is the number of real roots of /. But we can do better by applying LemmalU 



Lemma 7. Let r = maxjp, Tj,Tg}. Assume that we have isolated the r real roots of f and we know the signs 
of f over the isolating endpoints. Then, we can compute the sign of g over all r roots of f in Ob{p^(1t). 

Proof. Let Sj be the bitsize of the j'-th endpoint, where < j < r. The evaluation of SR(/, g) over 
this endpoint, by CoroUary |3l costs OsipciT + min{p, g}^Sj). To bound the overaU cost, we sum over 
ah isolating points. The first summand is OBip^qr). By Proposition [HI the second summand becomes 
OsCmiiijp, '?}^(p^ + pTf)) and is dominated. D 

3 Multivariate polynomials 

In this section, we extend the results of the previous section to multivariate polynomials, using binary 
segmentation ^. Let f,g € {Z[yi, . . . ,yk])[x] withdg^(/) = p > q = dg,^{g), dgy^{f) < di and dgy.{g) < dt. 

Let d — 11^=1 ^i ^^"^ ^{f) T^ {9) ^ '''• The y^-degree of every polynomial in SR(/, g) is bounded by 
dg^. (res(/, g)) < {p + q)di. Thus, the homomorphism ij} : Z[?/i, . . . , y^] -^ ^[y], where 

2/1^ y,y2^ 2/(p+«)'^S...,yfcK^ y(p+9)^-^rfi-rf.-i^ 

allows us to decode res('ip{f),ip{g)) = -(/'(res(/, (/)) and obtain res(/, g). The same holds for every polyno- 
mial in SR(/, y). Notice that '4'{f),'4'{g) G (^[y])[^] have y— degree less or equal to {p + q)''~^d since, in the 
worst case, f or g contains a monomial of the form y^'^ y^^ ■ ■ ■ y^J' ■ Thus, dg (res(7/'(/), tp(g))) < {p + q^d. 

Proposition 8. \28t l We can compute SRQ(/, g), any polynomial in SR(/, g), and res(/, g) w.r.t. x in 
dB{q{p + qf+^dT). 

Lemma 9. We can compute SR(/, g) in O B{q{p + q)^^"^ dr) . 

Proof. Every polynomial in SR(/, g) has coefficients of magnitude bounded 2'^*^^"'"'')'^, for a suitable constant 
c, assuming r > lg(d). Consider the map x '■ "^[y] '~^ ^i where y i-t- 2^'^^^'^''^'^^ ^ and let (/) = -0 ° X ■ 
Z[2/i, y2 . . . , yfc] ^ Z. Then £ (</.(/)) , C {<t>{g)) <c{p + qf dr. Now apply Proposition [B 

In order to complete the computation we should recover the result from the computed sequence, that is 
to apply the inverse image of </>. The cost of this computation (almost linear w.r.t. the output) is dominated; 
which is always the case. D 

Theorem 10. We can evaluate SIl{f,g) at x — a where a £ Q U {cxo} and C{a) = a, in Osiqip + 
q)'^+^dmax{r, cr}). 

Proof. First we compute SRQ(/,5) in OB{q{p+q)'^~^^dT) (Proposition|Sl) , and then we evaluate the sequence 
over a, using binary segmentation. For the latter we need to bound the bitsize of the resulting polynomials. 

The polynomials in SR(/, g) have total degree in j/i, . . . , y^ bounded by (p + q) J2i=i '^i ^^^'^ coefficient 
bitsize bounded by {p + q)T. With respect to x, the polynomials in SR(/, 5) have degrees in 0{p), so 
substitution x = a yields values of size 0{p<t). After the evaluation we obtain polynomials in Z[j/i, . . . , yk\ 
with bitsize bounded by max{(p + q)T,pa} < {p + q) max{r, a}. 

Consider the map x ■ '^[y] ^^ ^i where y 1-^ 2''^(p^'^> "'^^^^'^■'^'\ for a suitable constant c. Apply the map 
(j) = tp o X to f,g. Now, C {4>{f)) , C {4>{g)) < cd{p + q)'' maxjr, a}. By Proposition [21 the evaluation costs 
dB{q{p + q)''+^dina.x{T,a}). D 

We obtain the following, for bivariate /,5 G (Z[y])[a::], such that dg^(/) = p, dg^{g) = g, dgy{f),dgy{g) < 
d. 

Corollary 11. We compute S'R{f,g) in OB{pq{p + q)^dT). For any polynomial SRj(/,g) in SR(/,g), 
dg,(SR,(/,g)) = 0(max{p,g}),, dg,y(SR,(/,g)) = 0(max{p, g}d), and a/so £ (SR,(/,5)) - 0(max{p,q}T). 

Corollary 12. We compute SRCl{f,g), any polynomial in SIi{ f,g), and res{f, g) m ©^(pgmaxlp, gjdr). 

Corollary 13. We can compute SR(/, g ; a), where a £ QU{cxd} and C (a) = cr, in Os (pg inax{p, g}(imax{r, a}). 
For the polynomials SRj(/, (7; a) £ ^[y], except for f,g, it holds dg (SRj(/, g; a)) — 0{{p + q)d) and 
^ (SRj (/, g ; a)) = 0{ma.x{p, g}r + min{p, q}a) . 



Algorithm 1: SIGN_At(F, a,^) 


Input: F eZ[x,y],a^ {A,[ai 


32]),/? -(S, 


[bi,b2]) 


Output: sigii(i^(a,^)) 








1 compute SRQ^. (A, F) 








2 Li^SR,(A,F;ai), ^1^0 








3 foreach / e Li do Vi ^ ADD(yi, 


sign_at(/,/3)) 


4 L2^SR,iA,F;a2), V2 ^ 9 








5 foreach f & L2 do V2 ^ ADD(y2, 


sign_at(/,^)) 


6 RETURN (vAR(Vi) - VAR(V2)) • 


sign 


(A'ia)) 





We now reduce tlie computation of the sign of F e Z[x,y] over (a,/3) g R^, to that over several points 
in Q2. Let dg^(F) = dgy{F) = m, C{F) = a and a ^ (A, [31,82]), /3 = (5,[bi,b2]), where A,B e Z[X], 
dg{A) = dg(_B) = n2j_C {A) = £ (_B) = a. We assume rii < 77.2, which is relevant below. The pseudo-code is 



in Algorithm [U see [3l|, and generalizes the univariate case, e.g. 15|,|38|. For A, resp. B, wc assume that 
wc know their values on ai, 32, rcsp. bi, b2. 

Theorem 14. We compute the sign of polynomial F{x,y) over a, P in OB{n\n^(T). 

Proof. First, we compute SRQ^(A,F), in OBinln^cr) (Corollary [T2|). so as to evaluate SR(A, F) on the 
endpoints of a. 

We compute SR(A, F;3i). The first polynomial in the sequence is A and notice that we already know 
its value on ai. This computation costs OB{n\ rij cr) by CoroUarv 1131 with q = ni, p = 77.2, d ~ ni, t ~ a, 
and a = n2cr, where the latter corresponds to the bitsizc of the endpoints. After the evaluation we obtain a 
list Fi, which contains 0{ni) polynomials, say / G Z[y], such that dg(/) = 0(711712). To bound the bitsize, 
notice that the polynomials in SR(/, g) are of degrees 0{ni) w.r.t. x and of bitsize 0{n2(y). After we 
evaluate on Si, £(/) = 0{nin2(j). 

For each / € Li we compute its sign over /? and count the sign variations. We could apply directly 
Corollary [5l but we can do better. If dg(/) > 712 then SR(F, /) = [B, /, — F, g = — prem (/, — F) ,...). 
We start the evaluations at g: it is computed in OBinin^cr) (Proposition [T|) , dg(g) = 0(712) and C{g) = 
0(7ii7i2cr). Thus, we evaluate SR(— F,g;3i) in 05(7717420'), by Corollary [Sj with p = q ~ 7t2, t/ = a, 
T = 77i7i2cr. If dg(/) < 7*2 the Complexity is dominated. Since we perform 0(7ii) such evaluations, all of 
them cost OB{n\n2<j)- 

We repeat for the other endpoint of a, subtract the sign variations, and multiply by sign(A'(a)), which 
is known from the process that isolated a. If the last sign in the two sequences is alternating, then 
sign(F(a,/?)) = 0. D 

4 Bivariate real solving 

Let F,G & Z[x,y], dg(/) = dg{g) = 7i and C{F) = C{G) = a. We assume relatively prime polynomials 
for simplicity but this hypothesis is not restrictive because it can be verified and, if it does not hold, it 
can be imposed within the same asymptotic complexity. We study the problem of real solving the system 
F = G = 0. The main idea is to project the roots on their x- and y-coordinates. The difference between the 
algorithms is the way they match coordinates. 

4.1 The GRID algorithm 



Algorithm grid is straightforward, see also 1J,|36|. The pseudo-code is in Algorithm[2] We compute the x- 
and j/-coordinatcs of the real solutions by solving resultants reSa;(F, G), reSj,(F, G). We match them using 
the algorithm SIGN_AT (Theorem [H)) by testing all rectangles in this grid. 

To the best of our knowledge, this is the first time that the algorithm's complexity is studied. Its simplicity 
makes it attractive; however, SIGN_AT (Algorithm [1]) is very costly. The algorithm requires no genericity 
assumption on the input; we study a generic shear that brings the system to generic position in order to 



Algorithm 2: GKii){F,G) 



Input: F,G eZ[x,y] 

Output: The real solutions oi F = G = 

1 R,,^reSyiF,G) 

2 Lx.Mx ^— SOLVE(i?2;) 

3 Ry <— resi,(F, G) 

4 Ly,My ^ SOLVE(i?y) 

5 0^0 

6 foreach a e Lj; do 



foreach /3 G iy do 

|_ if SiGN_AT(i^, a,/3) = A SIGN_At(G, a,/3) = then Q ^ ABB{Q,{a, P}) 



9 RETURN Q 



compute the multiplicities within the same complexity bound. The algorithm allows the use of heuristics, 
such as bounding the number of roots, e.g. Mixed Volume, or counting the roots with given abscissa by 
Lemma [?!] 

Theorem 15. Isolating all real roots of system F = G = using grid has complexity Osin^^ + n^^a), 
provided a = 0{n^); or in OsiN^'^), where N — max{n,cr}. 

Proof. We compute resultant R^ of _F, G w.r.t. y (linc[T]in Algorithmic]). The complexity is OBin^a), using 
Corollary [T^ with p = q = d = n and t = a. Notice that dg{Rx) = 0{n^), C (Rx) = 0{na-). We isolate its 
real roots in Osin^^ + n^'^a'^) (Proposition 2]) and store them in L^- This complexity shall be dominated. 
We do the same for the y axis (lines [3] and |4] in Algorithm [2]) and store the roots in Ly. 

The representation of the algebraic numbers contains the square- free part of Rx or Ry, which has the 
bitsize ©(n^ + na) [3|, ll5[. The isolating intervals have endpoints of bitsize 0{n'^ + n^a). Let r^, ry be 
the number of real roots of the corresponding resultant, both in 0{n'^). For every pair of algebraic numbers 
from Lx and Ly, we test whether F, G vanish using SIGN_AT (Theorem [Til and Algorithm [T]) . Each test costs 
Osin^'^ + n^a) and we perform rx ry = 0(n'*) of them. D 

We now examine the multiplicity of a root {a, P) of the system. Refer to [J, Section IL6] for its definition 
as the exponent of factor (/3a; — ay) in the resultant of the (homogenized) polynomials, under certain as- 



sumptions. Previous work includes [16|, |32, |37[ . Our algorithm reduces to bivariate sign determination and 
does not require bivariate factorization. The sum of multiplicities of all roots (a,/3j) equals the multiplicity 
of a: = a in the respective resultant. We apply a shear transform so as to ensure that different roots project 
to different points on the x-axis. 

4.1.1 Deterministic shear 

We determine an adequate (horizontal) shear such that 

Rt{x) = resy {F{x + ty, y), G{x + ty, y)) , (1) 

has simple roots corresponding to the projections of the common roots of the system Fix, y) = G{x, y) = 0, 
when t !—>■ to € 'E, and the degree of the polynomials remains the same. Notice that this shear does not affect 
inherently multiple roots, which exist independently of the reference frame. Rred G (^M)[^] is the squarefree 
part of the resultant, as an element of UFD (Z[t])[a:], and its discriminant, with respect to x, is A G Z[t]. 
Then to must be such that A(io) 7^ 0. 

Lemma 16. Computing to G Z, such that the corresponding shear is sufficiently generic, has complexity 
(5B(nio + nV). 



Proof. Suppose to is such that the degree does not change. It suffices to find, among n'^ integer numbers, 
one that does not make A vanish; note that all candidate values are of bitsize O{\ogn). 

We perform the substitution {x, y) i— >■ (x+ty, y) to F and G and compute the resultant w.r.t. y in ©^(n^cr), 
which lies in Z[t, a:], of degree O^n?) and bitsize 0{da) (Proposition [8]). We consider this polynomial as 
univariate in x and compute its square-free part, and then the discriminant of its square-free part. Both 
operations cost Osin^^ 4- n^a) and the discriminant is a polynomial in Ti\t] of degree 0{n'^) and bitsize 
d{d'^ + (Pa) (Corollary [nD. 

We can evaluate the discriminant over all the first n* positive integers, in Osin^ + »T-^cr), using the 



multipoint evaluation algorithm, see e.g. |34|. Among these integers, there is at least one that is not a root 



of the discriminant. D 

The idea here is to use explicit candidate values of to right from the start. In practice, the above 
complexity becomes ©^(ri^cr), because a constant number of tries or a random value will typically suffice. 
For an alternative approach, see 17 1, and [3|. It is straightforward to compute the multiplicities of the 



sheared system. Then, we need to match the latter with the roots of the original system, which is nontrivial 
in practice. 

Theorem 17. Consider the setting of Theorem \15l Having isolated all real roots of F ^ G = 0, it is possible 
to determine their multiplicities in Osin^^ + n^^a + n^^a'^). 

Proof. By the previous lemma, t G Z is determined, with C (t) = 0{\ogn), in Osin^'^+n^o-). Using this value, 
we isolate all the real roots of Rt{x), defined in ([T]), and determine their multiplicities in Osin^^ + n^'^a'^) 
(Proposition HI). Let pj ~ {Rt{x), [rj,rj]) be the real roots, for j = 0, . . . , r — 1. 

By assumption, we have already isolated the roots of the system, denoted by (ai,/3i) G [ai,a^] x [6^,6^], 
where a^, a^, &i, &^ G Q for i = 0, . . . , r— 1. It remains to match each pair {on, pi) to a unique pj by determining 
function </) : {0, . . . , r — 1} — > {0, . . . , r — 1}, such that 0(i) = j iff {pj,f3i) G K^; is a root of the sheared 
system and a^ = pj + tj3i . 

Let [ci, c';] = [oi, a^] — t[&i, h'^ G Q^. These intervals may be overlapping. Since the endpoints have bitsize 
0{n'^ + n^a), the intervals [ci,c^] are sorted in Osin^ + n^a). The same complexity bounds the operation 
of merging this interval list with the list of intervals [rj,r'A. If there exist more than one [ci, c'^ overlapping 
with some [rj, r'], some subdivision steps are required so that the intervals reach the bitsize of Sj, where 2"^ 
bounds the separation distance associated to the j-th root. By Proposition [6l J2i ^i ~ C("-^ + n^a). 

Our analysis resembles that of [15| for proving Proposition H) The total number of steps is 0{^^ Si) ~ 
0{n'^ -\- n'^cr), each requiring an evaluation of R{x) over an endpoint of size < Sj. This evaluation costs 
Osin^Si), leading to an overall cost of Osin^ + rJ cr) per level of the tree of subdivisions. Hence, the overall 
complexity is bounded by Osin^^ + n^^a + n^"<T^). D 

4.2 The M_RUR algorithm 

M_RUR assumes that the polynomials are in Generic Position: different roots project to different cc-coordinates 
and leading coefficients w.r.t. y have no common real roots. 

Proposition 18. fj. \la] Let F,G be co-prime polynomials, in generic position. IfS'Rj{x,y) = sr j{x)y^ + 
srj,j^i{x)y^^^ + • • ■ -|- srjfi{x), and {a, /3) is a real solution of the system F = G = 0, then there exists k, 
such that sro(a) = • ■ • = srfe_i(Q) = 0, svk{a) ^ and l3 = -^^^j^^^^- 

This expresses the ordinate of a solution in a Rational Univariate Representation (RUR) of the abscissa. 
The RUR applies to multivariate algebraic systems [3|, |6|, 123 . |30[; by generalizing the primitive element 
method by Kroneckcr. Here we adapt it to small-dimensional systems. 



Our algorithm is similar to [16l [17[. However, their algorithm computes only a RUR using Proposition 
1181 so the representation of the ordinatcs remains implicit. Often, this representation is not sufficient (we 
can always compute the minimal polynomial of the roots, but this is highly inefficient). We modified the 



algorithm [IJ], so that the output includes isolating rectangles, hence the name modified- RUR (m_rur). 
The most important difference with [16| is that they represent algebraic numbers by Thom's encoding while 
we use isolating intervals, which were thought of having high theoretical complexity. 



Algorithm 3: m_rur {F, G) 


Input: F,G G Z[X, Y] in generic position 




Output: The real solutions of the system F = G = 




1 SR^SRj,(F,G) 




/* Projections and real solving with multiplicities 


*/ 


2 R.^^reSy{F,G) 




3 Px,M^ <- SOLVE(i?3;) 




4 Ry ^Tes^{F,G) 




5 Py,My ^ SOLVE{Ry) 




6 / ^ INTERMEDIATE_POINTS(/'y) 




/* Factorization of Rr^ according to sr 


*/ 


7 K ^ compute_k(SR, P^) 




8 Q^0 




/* Matching the solutions 


*/ 


9 foreach a £ P^ do 




10 


P <— FIND(a, K, Py,I) 




11 


_ Q^ADD(g,{«,/3}) 




12 RETURN Q 





The pseudo-code of M_RUR is in Algorithm [3l We project on the x and the y-axis; for each real solution 
on the X-axis we compute its ordinate using Proposition [181 First we compute the sequence SR(-F, G) w.r.t. 
y in OBirv'o) (Corollary [n]). 

Projection. This is similar to GRID. The complexity is dominated by real solving the resultants, i.e. 
Cb('i^^ + n^^ a^). Let cti, rcsp. /3j, be the real root coordinates. Wc compute rationals qj between the /3j's 
in ©^(n^cr), viz. iNTERMEDiATE_POiNTS(Fy); the q-j have aggregate bitsize 0{ii? a) (Lemma[6]): 

go < /3i < gi < /32 < • ■ • < fii-i < qe-i < l^e < qe, (2) 

where £ < 1r? . Every /3j corresponds to a unique a^. The multiplicity of a.i as a root of Rx is the multiplicity 
of a real solution of the system, that has it as abscissa. 

Sub-algorithm COMPUTE_k. In order to apply Proposition [TSl for every ai we must compute fc e N* 
such that the assumptions of the theorem are fulfilled; this is possible by genericity. We follow [iGJ, |25[ 
and define recursively polynomials V ^{x): Let $0(2^) = gcd(sro^x)?sr'(x)) ' *j(^) ^ gcd($j-i(a;),srj(x)), and 
Fj = -f-7;jY^, for J > 0. Now sri(x) e U\x\ is the principal subresultant coefhcient of SR; G (Z[x])[y], and 
<&o(a;) is the square-free part of Rx = sro(a;). By construction, $0(2;) = IljFjXa:^) and gcd(Fj,Fi) = 1, if 
J 7^ i. Hence every ai is a root of a unique Fj and the latter switches sign at the interval's cndpoints. Then, 
sro(a) = sri(a) = 0, . . . ,srj(a) = 0, srj-|_i(a) 7^ 0; thus k = j + 1. 

It holds that dg($o) = Oin^) and C ($0) = ^("^ +na). Moreover, J2j dg(Fj) = J^j ^j = C>{n^) and, by 
Mignotte's bound [2l[, C (Fj) — 0{n'^ + na). To compute the factorization $0(2;) = 111 ^ji^) ^^ ^ product 
of the srj(a;), we perform ©(n) gcd computations of polynomials of degree ©(n^) and bitsize 0{n^ + na). 
Each gcd computation costs Osin^ + n^ cr) (Proposition [T]) and thus the overall cost is OB{n^ + n^ <t). 

We compute the sign of the Fj over all the 0{n?) isolating endpoints of the a^, which have aggregate 
bitsize 0(n^ -|- n^ a) (LemmalB]) in OsiSju'^ + SjU^a + 5^(n* -I- n^cr)), using Horner's rule. Summing over all 
Sj, the complexity is Osin^ + rJcr). Thus the overall complexity is OBin^ + n^ a). 

Matching and algorithm FIND. The process takes a real root of Rx and computes the ordinate [3 
of the corresponding root of the system. For some real root a of Rx we represent the ordinate A{a) = 
— ^ ^'^s'r7(a) ~ XTaT- "^^^ generic position assumption guarantees that there is a unique /3j, in Py, such that 



Algorithm 4: G_RUR (F, G) 



Input: F,G gZ[X,Y] 

Output: The real solutions of the system F = G = 

/* Projections and real solving with multiplicities */ 

1 Ra: ^ reSy{F,G) 

2 P^,M^ ^solve{R^) 
s Ry ^ reSj.(F, G) 

4 Py,My <- SOLVE(i?y) 

/* / contains the rationals qi < q2 < ■ ■ ■ < q\j\ */ 

5 / ^ INTERMEDIATE_POINTS(Py) 

6 Q<-0 

7 foreach a & Px do 

8 F <— SQUAREFREEPART(F(a,y)) 

9 G <— SQUAREFREEPART(G(a,y)) 

10 H ^ GCd(F,G) e (Z[a])[j/] 

11 for j ^— 1 to |/| — 1 do 

12 if H{a, qj) ■ H{a, qj+i) < then 

/* Py[j] indicates the j-th element of Py */ 

13 \_Q^ADD{Q,{a,Py[j]}) 

14 RETURN Q 



I3j = A{a), where I < j < i- In order to compute j we use ^■. qj < A{a) = X(q) ^ 1^3 "^ 9i+i- Thus j can 
be computed by binary search in 0{\gi) = 0{\gn) comparisons of A{a) with the qj. This is equivalent to 
computing the sign of Bj{X) — Ai{X) — qj A2{X) over a by executing 0{\gn) times, SiGN_AT(i?j,a). 

Now, C{qj) = 0{n^ + nV) and dg(Ai) = dg(srfc,fc_i) = ©(n^), dgiA^) = dg(sr,.) = ©(n^), /:(Ai) = 
©(ncr), £(^2) = 0{n(j). Thusdg(Bj) = ^(ri.^) a.nAC{Bj) = 0{ri^+n^ a). We conclude that SIGN_At(Bj, a) 
and FIND have complexity 0B(n* + n^(T) ( Corollary [S|). As for the overall complexity of the loop (Lines CTTTj) 
the complexity is OB{n^^ + n^u), since it is executed 0{ii?') times. 

Theorem 19. We isolate all real roots of F = G = 0, if F, G are in generic position, by m_rur in 
OBin^'^ + ni°o-2-) . Q^ j„ dBiN^"^), where N = maxjri, a}. 

The generic position assumption is without loss of generality since we can always put the system in such 
position by applying a shear transform; see Section 14.1.11 and also [3|, |l6[ . The bitsize of polynomials of the 
(sheared) system becomes (D{n + a) [16[ and docs not change the bound of Theorem [T9l However, now is 
raised the problem of expressing the real roots in the original coordinate system (see the proof of Theorem 

4.3 The G_RUR algorithm 

In this section we present an algorithm that uses some ideas from m_rur but also relies on GCD computations 
of polynomials with coefficients in an extension field to achieve efficiency (hence the name G_RUr). The 
pseudo-code of G_RUR is in Algorithm |4l For GCD computations with polynomials with coefficients in an 
extension field we use the algorithm, and the MAPLE implementation, of van Hoeij and Monagan [33| . 

The first steps are similar to the previous algorithms: We project on the axes, we perform real solving 
and compute the intermediate points on the y-axis. The complexity is ©^(n^^ + n^^u'^). 

For each x-coordinate, say a, we compute the square- free part of F{a, y) and G(a, y), say F and G. The 
complexity is that of computing the gcd with the derivative. In [33[ the cost is OsirnMND + niN'^D'^ A- 
m?kD), where M is the bitsize of the largest coefficient, N is the degree of the largest polynomial, D is the 
degree of the extension, k is the degree of the gcd, and m is the number of primes needed. This bound does not 
assume fast multiplication algorithms, thus, under this assumption, it becomes OBimMND+mND+mkD). 
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In our case M = 0{a), N = 0{n), D = 0{n'^),k = 0{n), and m = 0{na). The cost is OBirfa'^) and 
since we repeat it 0{n^) times, the overall cost is ©^(n^cr^). Notice the bitsize of the result is Osin + o") 

i- 

Now for each a, we compute H = gcd(F, G). We have M = 0{n + <7), N = 0{n), D = ©(n^), k = 0{n), 

and m = 0{n^ +ncr), so the cost of each operation is Osin^ +7T."'cr^) and overall Osin^ +11^0-'^). The size of 
TO comes from Mignotte's bound |2l| . _ff is a square-free polynomial in (Z[a])[j/], of degree 0{n) and bitsize 
0{n^ + na), whose real roots correspond to the real solutions of the system with abscissa a. The crux of 
the method is that H changes sign only over the intervals that contain its real roots. To check these signs, 
it suffices to substitute y in _ff by the intermediate points, thus obtaining a polynomial in Z[a], of degree 
0{n) and bitsize 0{n^ + na + nsj), where Sj is the bitsize of the j-th intermediate point. 

Now, we consider this polynomial in Wi[x] and evaluate it over a. Using Corollary [5] with p = n^, 
Tf ~ n^ + na, q = n, and Tg = n'^ + na + nsj, this costs Osin^ + n^cr + n'^Sj). Summing over 0{n'^) points 
and using Lemma [SI we obtain Osin^ + n^cr). Thus, the overall complexity is 05(71^" + n^a). 

Theorem 20. We can isolate the real roots of the system F = G = 0, using G_RUR in Osin^^ + n^'^cr^); or 
OsiN^^), where N ~ maxjri, a}. 

5 Applications 

5.1 Real root counting 

Let F G Z[a;, y], such that dg^{F) — dgy{F) = 711 and jC {F) = a. Let a,j3e ^aig, such that a = (A, [ai, 82]) 
and 13 = {B, [bi, bs]), where dg(A), dg(B) = na, £ (A) ,C{B) <t and c e Q, such that C (c) = A. Moreover, 
assume that n\ = 0(77-2), as is the case in applications. We want to count the number of real roots of 
F = F{a,y) € {Z{a))[y] in (—00, +00), in (c, +00) and in {(3, +00). We may assume that the leading 
coefficient of F is nonzero. This is w.l.o.g. since we can easily check it, and/or we can use the good 
specialization properties of the subresultants [16|, |l8|, |20[ . 

Using Sturm's theorem, e.g. [l,!!!], the number of real roots of f' is VAR(SR(F, Fy-, -oo))-VAR(SR(^, Fy-, +00)). 
Hence, we have to compute the sequence STt{F,Fy) w.r.t. y, and evaluate it on ±00 or, equivalently, to 
compute the signs of the principal subresultant coefBcients, which lie in Z(a). This procedure is equivalent, 
due to the good specialization properties of subresultants |3|, |l8| , to computing the principal subresultant 
coefficients of S'R.{F,Fy), which are polynomials in Z[a;], and to evaluate them over a. In other words, the 
good specialization properties assure us that we can compute a nominal sequence by considering the bivariate 
polynomials, and then perform the substitution x = a. 

The sequence sr of the principal subresultant coefficients can be computed in OB{nfa), using Corollary 
[T^lwith p = q = d = ni, and t = a. Now, sr contains 0{ni) polynomials in Z[.t], each of degree 0{n1) and 
bitsize 0{nia). We compute the sign of each one evaluated over a in 

Ob(71]^712 maxJT, nia} + 712 uiin{ni, 712} r) 
using Corollary [S] with p = 712, q = nf, Tf = t, and Tg = nia. This proves the following: 
Lemma 21. We count the number of real roots of F ^ F{a,y) in OB{nfn2a + 77^7721"). 

In order to count the real roots of F in (/3, +00), we use again Sturm's theorem. The complexity of the 
computation is dominated by the cost of computing VAR{STl{F,Fy; (3)), which is equivalent to computing 
SR{F, Fy) w.r.t. to 7/, which contains bivariate polynomials, and to compute their signs over (a, /3). The 
cost of computing SR(i^, Fy) is OBinfa) using Corollary [TT] with p = q ^ d = ni, and t = cr. The sequence 
contains 0(771) polynomials in Z[a;,7/] of degrees 0(771) and 0(77i), w.r.t. x and y respectively, and bitsize 
0{nia). We compute the sign of each over (a,/3) in Os(77f772 max{7iicr, r}) (Theorem [H)) . This proves the 
following: 

Lemma 22. We count the number of real roots of F in {l3, +00) in 03(77^772 maxjnicr, r}). 
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By a more involved analysis, taking into account the difference in the degrees of the bivariate polynomials, 
we can gain a factor. We omit it for reasons of simplicity. Finally, in order to count the real roots of F in 
(c, +00), it suffices to evaluate the sequence S'R{F^ Fy) w.r.t. y on c, thus obtaining polynomials in Z[a;], 
and compute their signs over a. ^ 

The cost of the evaluation SR(F,i^j,;c) is C'b(»t-i uiaxjcr, A}), using CoroUarv fT51 with p = q = d = rii, 
T ~ a and a = X. The evaluated sequence contains 0{ni) polynomials in Z[a;], of degree 0{n1) and bitsize 
0{ni maxjcr, A}). The sign of each one evaluated over a can be computed in 

OB{nin2 maxJT, nicr, niA} + n^n2T), 
using Corollary [S] with p = 71,2, q^ n\, Tf = t and Tg = rti maxja. A}. This leads to the following: 
Lemma 23. We count the number of real roots of F in (c, +00) in OB(n|n2 max{niT, cr, A}). 

5.2 Simultaneous inequalities in two variables 

Let P,Q, Ai, . . . , Ai-^ , Bi, . . . , Bi^ , Ci , . . . , C^^ £ Z[X, Y] , such that their total degrees are bounded by n and 
their bitsize by a. We wish to compute (a, (3) G R^, such that P{a, /3) = Q{a, (3) = and also Ai{a, (i) > 0, 
Bj{a,P) < andCfc(a,^) = 0, where 1 < i < ii,l < j < £2,1 < k < £3. Let e = ii + £2 + is- 

Corollary 24. There is an algorithm that solves the problem of i simultaneous inequalities of degree < n 
and bitsize < a, in Osiin^^ + £n^^a + n^^a^). 

Proof. Initially we compute the isolating interval representation of the real roots of P = Q = in Osin^'^ + 
n^^a^), using GRUR_SOLVE. There are 0{n'^) real solutions, which are represented in isolating interval 
representation, with polynomials of degrees 0{n^) and bitsize 0{n'^ + no). 

For each real solution, say (a, /3), for each polynomial Ai, Bj, Ck we compute the signs of sign(Ai(a, /3)), 
sign {Bi{a, /?)) and sign {Ci{a, f3)). Each sign evaluation costs ©^(n^^ + n^cr), using Theorem [HI with rii = n, 
n2 = n^ and a = n^ + na. In the worst case we need n^ of them, hence, the cost for all sign evaluations is 
OB(£ni2 + £n"CT). D 

5.3 The complexity of topology 

In this section we consider the problem of computing the topology of a real plane algebraic curve, and improve 



upon its asymptotic complexity. The reader may refer to e.g. [3|, [ly, |25| for the details of the algorithm. 

We consider the curve in generic position fSection l4.1.ip . defined by F S Z[x, y], such that dg(/) = n and 
C{F) = a. We compute the critical points of the curve, i.e. solve F = _F!y = in Osin^^ + n^'^a'^), where 
Fy is the derivative of F w.r.t. y. Next, we compute the intermediate points on the x axis, in ©^(n* + n'^a) 
(Lemma [5]). For each intermediate point, say g^, we need to compute the number of branches of the curve 
that cross the vertical line x = q.j. This is equivalent to computing the number of real solutions of the 
polynomial F(qj, y) £ Z[y], which has degree d and bitsize 0{nC {q-j))- For this we use Sturm's theorem and 
Theorem [21 and the cost is Osin'^^ {<lj))- For all q^-'s the cost is Osin^ + n^cr). 

For each critical point, say (a, /3) we need to compute the number of branches of the curve that cross the 
vertical line x = a, and the number of them that are above y = /3. The first task corresponds to computing 
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the number of real roots of F{a, y), by application of Lemma \2T\ in Osin^ + "-^o"), where ni = n, 712 = n 
and T = n^ + na. Since there are 0{n'^) critical values, the overall cost of the step is Osin^^ + n^^a). 

Finally, we compute the number of branches that cross the line x = a and are above y = /3 in Osin^^ + 
n^'^a), by Lemma [22] Since there are 0{n'^) critical points, the complexity is Osin^^ +n^'^a). It remains to 
connect the critical points according to the information that we have for the branches. The complexity of 
this step is dominated. It now follows that the complexity of the algorithm is Osin^^ + n^^a + n^^tr^), or 
Ob{N'^^), which is worse by a factor than [3[. 

We improve the complexity of the last step since m_rur computes the RUR representation of the ordi- 
natcs. Thus, instead of performing bivariate sign evaluations in order to compute the number of branches 
above y = /3, we can substitute the RUR representation of /3 and perform univariate sign evaluations. This 
corresponds to computing the sign of 0{n^) polynomials of degree 0{n?) and bitsize 0{n'^ + n^a), over all 
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the a's [16|. Using Lenima[7]for each polynomial the cost is Osin^^ + n^cr), and since there are Os(n^) of 
them, the total cost is Osin^^ + ?i^^cr). 

Theorem 25. We compute the topology of a real plane algebraic curve, defined by a polynomial of degree n 
and hitsize a, in Osiri^^ + n^^cr + n^^a^), or Ob{N^^), where N = max{n,cr}. 

Thus the overall complexity of the algorithm improves the previously known bound by a factor of N'^. 
We assumed generic position, since we can apply a shear to achieve this, see Section [JT] 

6 Implementation and Experiments 

This section describes our open source maple implementation [^ and illustrates its capabilities through 
comparative experiments. Refer to [10[ for its usage and further details. Our design is object oriented and 
uses generic programming in view of transferring the implementation to C++ in the future. 

We provide algorithms for signed polynomial remainder sequences, real solving of univariate polynomials 
using Sturm's algorithm, computations with one and two real algebraic numbers, such as sign evaluation 
and comparison and, of course, solving bivariate systems. 

6.1 Our solvers 

The performance of all algorithms is averaged over 10 executions on a maple 9.5 console using a 2GHz 
AMD64@3K+ processor with 1GB RAM. The polynomial systems tested are given in [10[: systems Ri,Mi, Di 
arc from [l4|, the Ci are from [17|, and Wi,i = 1,...,4, follow from CiS after swapping x,y. The latter 



arc of the form / = ^ = 0. For gcd computations in a (single) extension field, the package of [33| is used. 



dy 

The optimal algorithms for computing and evaluating polynomial remainder sequences have not yet been 
implemented. 

Our main results are reported in Table [1] G_RUR is the solver of choice since it is faster than grid and 
M_RUR in 17 out of the 18 instances. However, this may not hold when the extension field is of high degree. 
G_RUR yields solutions in < 1 sec, apart from C5. For total degree < 8, G_RUR requires < 0.4 sec. On 
average, G_RUR is 7-11 times faster than grid, and about 38 times faster than m_rur. The inefficiency of 
M_RUR is due to the fact that it solves sheared systems which are dense and of increased bitsizc; it also 
computes multiplicities. Finally, GRID reaches a stack limit with the default maple stack size (8, 192 KB) 
when solving C5. Even when we multiplied stack size by 10, grid did not terminate within 20 min. 

Whenever we refer to the speedup we imply the fraction of runtimes. G_RUR can be up to 21.58 times 
faster than grid with an average speedup of around 7.27 among the input systems (excluding C5). With 
respect to m_rur, G_RUR can be up to 275.74 times faster, with an average speedup of 38.01. 

Filtering has been used. For this, two instances of isolating intervals are stored; one for filtering, another 
for exact computation. Probably, the most significant filtering technique is interval arithmetic. When 
computing the sign of a polynomial evaluated at a real algebraic number, the first attempt is via interval 
arithmetic, applied along with [l|. When this fails, and one wants to compare algebraic numbers or perform 
univariate SIGN_AT, then the gcd of two polynomials is computed. 

Filtering helps most with m_rur, especially when we compute multiplicities. With this solver, one more 
filter is used: the intervals of candidate a;-solutions are refined by [l| so as to help the interval arithmetic 
filters inside find. If the above fails, we switch to exact computation via Sturm sequences, using the initial 
endpoints since they have smaller bitsize. In grid's case, filtering provided an average speedup of 1.51, 
where C5 has been excluded. With G_RUR, we have on average a speedup of 1.08. This is expected since 
G_RUR relies heavily on gcd's and factoring. 

Figure fTal shows the runtime breakdown corresponding to the various stages of each algorithm: Pro j ections 
shows the time for computing resultants, Univ. Solving for solving them, and Sorting for sorting solutions. 
In grid's and m_rur's case, biv. solving corresponds to matching. In G_RUR's case, matching is divided 
between rational biv and Raig biv; the first refers to when at least one of the co-ordinates is ratio- 
nal. Inter. points refers to computing intermediate points between resultant roots along the y-axis. StHa 
seq refers to computing the StHa sequence. Filter x-cand shows the time for additional filtering. Compute 
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phase of the 
algorithm 


interval 


median 


mean 


std 
dev 


mill 


max 


D 

5 
o 


projections 


00.00 


00.53 


00.04 


00.08 


00.13 


univ. solving 


02.05 


99.75 


07.08 


26.77 


35.88 


biv. solving 


O0.f9 


97.93 


96.18 


73.03 


36.04 


sorting 


00.00 


01.13 


00.06 


00.12 


00.26 




projection 


00.00 


00.75 


00.06 


00.14 


00.23 


univ. solving 


00. f8 


91.37 


15.55 


17.47 


20.79 


StHa seq. 


00.08 


38.23 


01.17 


05.80 


09.91 


inter, points 


00.00 


03.23 


00.09 


00.32 


00.75 
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5 
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9 


10 


2 
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389 
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Di 
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5 
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6 
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D2 


2 


2 


4 
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126 


Ci 
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6 


6 


1,702 
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C2 
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3 


6 
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99 


C,i 


8 


7 


13 


669 


1,815 


152 


Ci 


8 


7 


17 


7,492 
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474 
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16 


15 


17 
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6,367 


Wi 
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6 


9 
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W-i 


4 


3 


5 
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W'i 
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7 


13 


1,769 


2,333 


230 
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7 


17 


5,783 


77,207 
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(a) Statistics on SLV's sub-algorithms. 



(b) Performance of our solvers when computing mul- 
tiplicities. 



Figure 1: Statistics 



K reflects the time for sub-algorithm COMPUTE-k. In a nutshell, grid spends more than 73% of its time 
in matching. Recall that this percent includes the application of filters and does not take into account C5 . 
M_RUR spends 45-50% of its time in matching and 24-27% in filtering. G_RUR spends 55-80% of its time in 
matching, including gcd computations in an extension field. 

In order to compute multiplicities, the initial systems were sheared whenever it was necessary, based on 
the algorithm presented in Section [4.1.11 Overall results are shown in Figure [Tbl grid's high complexity 
starts to become apparent. Overall, G_RUR is fastest and terminates within < 1 sec. It can be up to 15.81 
times faster than grid with an average speedup of around 5.26. With respect to m_rur, this time G_RUR 
can be up to 170.15 times faster, with an average speedup of around 18.77 among all input polynomial 
systems. m_rur can be up to 6.23 times faster than GRID, yielding an average speedup of 1.71. A detailed 
table in [lO[ gives us the runtime decomposition of each algorithm in its major subroutines. Results are 
similar to Section l6Jl except that G_RUR spends 68-80% of its time in matching, including gcd's. In absence 
of excessive factoring G_RUR spends significantly more time in bivariate solving. 

6.2 Other software 

fgb/rs □ [30|] performs exact real solving using Grobner bases and RUR, through its MAPLE interface; 
additional tuning might offer 20-30% efficiency increase. Three SYNAPS |3 solvers have been tested: STURM 
is a naive implementation of grid IJ]; SUBDIV implements [23[, using the Bernstein basis and double 
arithmetic. It needs an initial box and [—10,10] x [—10,10] was used, newmac [2J] is a general purpose 
solver based on eigenvectors using lapack, which computes all complex solutions. 

MAPLE implementations: insulate implements 137| for computing the topology of real algebraic curves, 
and TOP implements |17| . Both packages were kindly provided by their authors. We tried to modify the 
packages so as to stop as soon as they compute the real solutions of the corresponding bivariate system. It 
was not easy to modify insulate and top to deal with general systems, so they were not executed on the 
first data set. top has a parameter that sets the initial precision (decimal digits). There is no easy way for 
choosing a good value. Hence, recorded its performance on initial values of 60 and 500 digits. 

Experiments are not considered as competition, but as a crucial step for improving existing software. 
It is very difficult to compare different packages, since in most cases they are made for different needs. In 



"http : //www-spaces . Iip6 . f r/index . html 

^http : //www-sop . Inria . f r/galaad/logiciels/synaps/ 
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Table 1: Performance of our solvers and other tested software. 
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G_RUR 


STURM 


SUBDIV 


NEWMAC 


60 


500 
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- 


- 
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1 


66 
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36 
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- 


- 
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1 


1 
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87 


72 
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2 


- 


- 


- 
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4 


2 


3 


4 


5 


4 


24 


1 


289* 


2 


- 


- 


- 


A/3 


6 


3 


5 


803 


782 


110 


30 


230 


5,058* 


7 


- 


- 


- 


A/4 


9 


10 


2 


218 


389 


210 


158 


90 


3* 


447 


- 


- 


- 


Di 


4 


5 


1 


6 


12 


6 


28 


2 


5 


8 


- 


- 


- 


D2 


2 


2 


4 


667 


147 


128 


26 


21 


1* 


2 


- 


- 


- 


Ci 


7 


6 


6 


1,896 


954 


222 


93 


479 


170,265* 


39 


524 


409 


1,367 


C2 


4 


3 


6 


177 


234 


18 


27 


12 


23* 


4 


28 


36 


115 


C3 


8 


7 


13 


580 


1,815 


75 


54 


23 


214* 


25 


327 


693 


2,829 


Ci 


8 


7 


17 


5,903 


80,650 


370 


138 


3,495 


217* 


190* 


1,589 


1,624 


6,435 


C5 


16 


15 


17 


>20' 


60, 832 


3,877 


4,044 


>20' 


6,345* 


346* 


179, 182 


91,993 


180,917 


w. 


7 


6 


9 


2,293 


2,115 


247 


92 


954 


55,040* 


39 


517 


419 


1,350 


W2 


4 


3 


5 


367 


283 


114 


29 


20 


224* 


3 


27 


20 


60 


W3 


8 


7 


13 


518 


2,333 


24 


56 


32 


285* 


25 


309 


525 


1,588 


W4 


8 


7 


17 


5,410 


77, 207 


280 


148 


4,086 


280* 


207* 


1,579 


1,458 


4,830 



addition, accurate timing in maple is hard, since it is a general purpose package and a lot of overhead 
is added to its function calls. Lastly, the amount of experiments is not very large in order to draw safe 
conclusions. 

Overall performance results are shown on Table [1] In cases where the solvers failed to find the correct 
number of real solutions we indicate so with *. Note that in newmac's column an additional step is required 
to distinguish the real solutions among the complex ones. In the sequel we refer only to G_RUR, since it is 
our faster implementation. 

G_RUR is faster than fgb/rs in 8 out of the 18 instances, including C5. The speedup factor ranges from 
0.2 to 22 with an average of 2.62. 

As for the three solvers from SYNAPS, G_RUR is faster than STURM in 6 out of the 18 instances, but it 
behaves worse usually in systems that are solved in < 100 msecs. because STURM is implemented in C++. As 
the dimension of the polynomial systems increases, G_RUR outperforms STURM and the latter's lack of fast 
algorithms for computing resultants becomes more evident. Overall, an average speedup of 2.2 is achieved. 
Compared with SUBDIV, G_RUR is faster in half of the instances and similarly to the previous case is slower 
on systems solved in < 400 msecs. On average, G_RUR achieves a speedup of 62.92 which is the result of the 
problematic behavior of SUBDIV in Ci and Wi. If these systems are omitted, then the speedup is 8.93 on 
average, newmac is slower than G_RUR in Mi,Di and W^ and comparable in Ri and R^. This time the 
average speedup of our implementation is 0.53. There are cases where newmac may not compute some of 
the real solutions. 

Finally, concerning the other maple software, insulate is slower than G_RUR in all systems but W2, 
thus our solver achieves an average speedup of 8.85. Compared to top with 60, resp. 500, digits, G_RUR is 
faster in all systems but Vl^2, yielding an average speedup of 7.79, resp. 22.64. Moreover, as the dimension 
of the polynomial systems increases, it becomes more efficient. 
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